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Chapter  1 


INTRODUCTION 

The  continuing  search  for  lighter,  stronger,  more  economical 
structures,  particularly  in  the  aircraft  and  missile  industry,  has  led 
to  the  investigation  of  various  composite  materials  as  a possible  appli- 
cable type  of  construction.  Many  of  these  materials  are  orthotropic  and 
multi-layered  and  have  nonlinear  physical  properties.  Tests  on  actual 
structures  have  proven  that  some  of  these  materials  are  unlike  conven- 
tional structural  materials  in  many  respects.  Although  much  has  been 
learned  about  composite  material  behavior  in  the  past  few  years,  there 
are  still  many  areas  that  are  unknown  and  unpredictable.  These  charac- 
teristics pose  many  difficult  problems  for  the  designer  and  analyst. 

To  obtain  the  most  efficient  use  of  these  composites,  suitable  analysis 
and  design  techniques  must  be  developed. 

The  aircraft  and  missile  industries  have  many  applications  for 
shells  of  revolution  and  are  continually  searching  for  ways  to  decrease 
the  weight  while  maintaining  or  increasing  the  strength.  The  possibility 
of  achieving  this  by  use  of  composites  has  created  much  interest  in  their 
application.  The  need  for  techniques  to  analyze  these  structures  was  the 
prime  motivation  for  this  present  effort.  This  study  investigates 
an  orthotropic,  laminated  shell  of  revolution  with  shear  deformation. 

The  rapid  development  of  the  digital  computer  in  the  last  two 
decades  has  contributed  significantly  to  the  analyst's  ability  to  treat 
these  complex  problems.  Numerical  methods  that  would  have  been  impossible 
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to  employ  prior  to  the  computer  are  now  quite  popular  and  are  used 
extensively  throughout  the  stress  analysis  community.  Both  the  finite 
element  method  and  finite  difference  techniques  are  used  today  to  solve 
complex  structural  problems.  Each  of  these  has  features  which  make  its 
application  more  suitable  to  a particular  type  problem. 

The  finite  element  method  was  chosen  for  this  present  development 
for  several  reasons.  Because  of  the  flexibility  of  their  size  and  shape 
finite  elements  can  represent  a given  body,  however  complex  its  shape 
may  be,  quite  accurately.  Structures  with  holes  or  discontinuities  can 
be  treated  with  little  difficulty.  Problems  involving  variable  material 
properties  and  geometry  such  as  are  encountered  with  fiber  composites  do 
not  present  any  additional  difficulty.  Geometrical  and  material  non- 
linearity can  be  dealt  with  relatively  easily.  One  of  the  principal 
assets  of  the  finite  element  method  is  the  ease  with  which  boundary 
conditions  can  be  represented. 

The  essential  feature  of  the  finite  element  method  is  that  the 
governing  differential  equations  of  equilibrium  of  the  shell  are  approxi 
mated  by  a set  of  algebraic  equations.  This  is  equivalent  to  substituti 
for  the  actual  structure,  an  assemblage  of  discrete  elements  inter- 
connected at  a finite  number  of  nodes.  The  element  stiffness  is  then 
evaluated  and  superimposed  to  obtain  the  total  stiffness  matrix  of  the 
entire  structure.  Finally,  the  nodal  force  equilibrium  equations  are 
solved  simultaneously  for  the  nodal  displacements  of  the  complete 
system  [ 1 -5 ] . 

In  shells  of  revolution,  the  structure  is  divided  into  a number 
of  short  frustums  which  are  connected  at  their  nodal  circle.  The 
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assemblage  is  made  through  equilibrium  and  compatibility  requirements 
at  the  nodal  circle.  Mayer  and  Harmon  [6]  employed  the  conical  frustum 
(singly  curved)  element  in  the  earlier  stage  of  analysis  of  shells  of 
revolution.  Popov  et  al.  [71  used  the  bending  displacements  due  to 
edge  loadings  from  the  exact  shell  theory  rather  than  the  usual  assumed 
displacement  functions.  Their  result  showed  no  significant  improvement 
over  the  simpler  assumed  functions.  Grafton  and  Strome  [8]  used  conical 
elements  in  a true  finite  element  technique.  The  results  are  very 
satisfactory  for  shells  with  a straight  generated  curve.  However,  there 
are  still  some  inaccurate  moments  due  to  the  approximation  of  a doubly 
curved  shell  by  a singly  curved  element.  This  is  mainly  caused  by  the 
discontinuity  of  slope  at  the  nodal  circle  of  the  substitute  structure. 

To  remedy  this  problem,  Jones  and  Strome  [9]  developed  a doubly  curved 
element  which  matched  both  the  location  and  slope  of  the  original  shell 
at  nodal  circles,  thus  avoiding  unwanted  discontinuities  of  slope  at 
these  locations.  Kho jasteh-Bakht  [10]  used  two  local  coordinate  systems, 
viz.,  curvilinear  and  rectilinear  coordinate  systems,  to  formulate  his 
element.  The  latter  proved  to  be  a better  approach  because  it  can 
treat  certain  constant  strain  states  which  the  former  was  unable  to 
acdommodate . 

Recently,  a f inite -element  technique  including  transverse  shear 
effect  has  been  attempted,  Clough  and  Felippa  [11]  have  described  a 
simple  shear  distortion  mechanism  which  can  be  implemented  by  expressing 
the  total  rotation  of  a cross  section  as  the  sum  of  the  rotation  on  the 
middle  surface  plus  a uniform  shear  strain  through  the  thickness. 

Klein  [12]  has  applied  the  matrix  displacement  finite  element 
approach  to  the  linear  elastic  analysis  of  shells  of  revolution  under 
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axisymmetr ic  loads.  The  shell  is  idealized  as  a series  of  conical 
frusta,  joined  at  nodal  circles.  The  external  forces  are  applied  at 
the  nodal  circles.  A comparison  with  these  solutions  is  made  in 
Chapter  7. 

Sharifi  [13]  used  an  incremental  formulation  for  the  nonlinear 
finite  element  analysis  of  sandwich  structures.  The  nonlinearities 
considered  were  due  to  large  displacements.  Included  in  the  analysis 
are  axisymmetr ic  shells  with  axisymmetric  loadings  and  boundary  condi- 
tions. Curved  elements  were  used  in  the  development, 

McNamara  [14]  investigated  nonlinear  dynamic  problems  by  using 
an  incremental  stiffness  finite  element  analysis.  Both  geometric  and 
material  non  1 inear it ies  were  considered. 

Becker  and  Brisbane  [15]  developed  the  equations  for  an  axisym- 
metric, orthotropic  shell  using  a straight  line  element.  Their  develop- 
ment did  not  include  shear  deformations.  Shear  deformations  are 
important  for  the  analysis  of  fiber  composites, 

Nickell  and  Sato  [16]  used  a curved  shell  element  to  analyze  an 
orthotropic,  layered  shell  of  revolution.  Shear  deformations  were  not 
included  in  their  analysis. 

In  this  present  effort,  the  finite  element  method  is  used  with 
a curved  shell  element  considering  a nonlinear  laminated,  orthotropic 
shell  of  revolution  and  transverse  shear  deformations. 

A polynomial  representation  of  the  meridional  curve  of  the  shell 
was  chosen  which  matches  the  displacements,  slope,  and  curvature  at  the 
nodal  points.  Nonlinear  terms  are  included  in  the  strain-displacement 
relationships.  The  stiffness  matrix  is  then  derived  using  these  non- 
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There  are  four  degrees  of  freedom  at  a node,  viz,,  two  translation, 
one  bending  rotation  and  one  transverse  shear  rotation.  The  field 
equations  similar  to  Reissner's  theory  of  thick  plates  [5]  were  used  as 
a guideline  for  formulating  the  shear  deformation  degree  of  freedom. 

The  procedure  employed  was  similar  to  that  of  Clough  and  Felippa  [11], 

The  classical  Kirchhof f -Love  assumption  for  normals  to  the  midsurface 
was  relaxed  in  favor  of  the  assumed  shear  deformation  mode. 

A computer  program  was  written  implementing  the  derived  equations. 
The  element  stiffness  matrix  was  formed  by  numerical  integration.  Much 
of  the  data  is  generated  internally  in  the  program.  The  program  is 
limited  to  ten  different  materials  and  50  nodes,  but  it  can  be  increased 
by  increasing  the  dimension  statement  accordingly.  The  program  is 
relatively  fast.  Most  of  the  example  problems  run  required  from  3 to  4 
seconds  execution  time  on  the  CDC  6600  computer. 


Chapter  2 


ELEMENT  GEOMETRY 

The  shell  to  be  considered  is  axisymmetr ic ; therefore,  it  is 
sufficient  to  define  only  the  shape  of  its  meridional  curve.  The  finite 
element  mechod  will  be  used  for  this  analysis.  The  element  is  shown  in 
Figure  1. 


Figure  1.  Curved  Shell  Element. 


The  element  is  curved  and  the  two  end  points  of  the  element  are 
denoted  by  I and  J.  For  a shell  whose  reference  surface  is  a surface  of 
revolution,  the  lines  of  principal  curvature  are  its  meridians  and 
parallel  circles.  The  principal  curvilinear  coordinates  of  the  reference 
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surface  are  the  angle  <P  between  the  normal  to  the  reference  surface  and 
the  axis  of  revolution  and  the  angle  d describing  the  position  of  points 
on  the  corresponding  parallel  circle.  Since  this  development  is  axisym- 
metric  both  in  geometry  and  load,  it  is  independent  of  d. 

A local  coordinate  system  for  each  curved  element  is  constructed 
between  two  adjacent  nodal  circles  by  drawing  chords  between  the  points. 
This  rectangular  Cartesian  system  which  is  normalized  by  the  chord 
length  £ is  denoted  by  x-y.  The  global  coordinates  are  represented  by 
r-z.  The  angles  shown  in  Figure  1 are  related  by  the  relation 

0 + i|f  + p = f (1) 


The  angle  4>  is  the  angle  between  the  normal  to  the  reference 
surface  and  the  axis  of  revolution.  The  angle  \|r  is  the  angle  between 
the  chord  of  the  element  and  the  z-axis.  The  angle  p is  defined  to  be 
the  angle  between  the  chord  line  (the  x-axis)  and  the  tangent  to  the 
curved  surface  at  any  point. 

From  Equation  (1), 

sin  P = cos  (<t>  + \|f)  = cos  <t>  cos  - sin  <t>  sin 

(2) 

cos  P = sin  (<t>  + \|c)  = sin  0 cos  \|r  + cos  <t>  sin  \|r 

To  approximate  the  meridional  curve,  the  following  substitute 
curve  is  assumed: 


(3) 


where  £ = chord  length  of  the  element. 


(0  < x < £) 


r 


Differentiating  Equation  (3)  with  respect  to  x 


£-*i  + 


2(a,  - a,)  3(a^  - a?)  2 4(a4  - a3>  3 4 

x + o x + o x " l X 


(4) 


df*  2(a2  - V . 6<‘l  - V , + 121 (a*  ~ ^ x2  . x3 

dx2  * 22  *3  « 


The  constants  a^  a.,,  a3>  and  a4  can  be  determined  by  evaluating 
Equations  (4)  at  the  end  points 


a^  = tan 


3r 

a2  = tan  P1  ' 2^  SeC  P1 
_3„  . 3 


(5) 


= sec3p2  + ~~  sec  ^ - 4 tan  ^ - 5 tan  pJ 


3 £ 3 

sec 


>4  = 2rT  S£C  P2  • 2Rl 


*P  + 3(tan  + tan  . 


where  the  subscripts  l and  2 on  p and  R reference  these  items  to  the 
I and  J nodes  respectively. 


The  following  geometrical  relations  are  given  with  respect  to  the 


element : 


Ar  = r2  - rL 


Az  = - z2 


^(Ar)2  + (Az)2 


1 Ar 

sin  V = — 


Az 

cos  l|l  = “ 


(6) 


T.: 
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After  the  substitute  curve  has  been  established,  all 
quantities  can  be  written  as  follows: 


tan  (3 


dx 

~ r ^ x sin  + y cos 


Since 


dr 

dx 


— sin  + tan  p cos  ^ 


— o d 

dS  = cos  P ^ 


i£  = <!£  dS  _ 1 

dx  dS  dx  ~ " R sec  P 


dy 

dJ  = tan  e 


therefore 


d y d . 

^2  = <tan  P) 


d£  _ i 


20  d6 
sec  6 "~ 
dx 


dx 


= - - sec  p 


The  quantity  dp/dS  is  negative  since  p is  decreasing  as  S is 
Therefore 


and 


d_y  1 3„ 

dx2  " « S“  13 


cos  4>  = sin  P cos  + cos  P sin  ^ 
sin  4>  = cos  p cos  \|r  - sin  P sin  \\i 


the  geometric 


(7) 


increasing . 


(8) 


(9) 


Chapter  3 


TRANSFORMATION  OF  COORDINATE  SYSTEMS 


The  displacement  vector  of  a material  point  on  the  midsurface 
in  the  local  principal  curvilinear  shell  coordinate  is  denoted  by 


where 


(f  } = [u  , w , 

c c c 


7c] 


u^  = the  displacement  along  the  meridian. 

= the  transverse  (normal)  displacement. 

< = the  rotation  about  a meridional  tangent. 

7 = shear  deformation, 

c 


(10) 


The  displacement  components  which  refer  to  the  local  rectilinear 
coordinates,  x-y,  are 


If  jT  = 

r 


X , 7 

r r 


(ID 


and  to  the  global  coordinates,  r-z,  are 

{f)T  = [u,  w,  X,  y)  . (12) 

The  transformation  between  these  components  can  be  seen  as 
fol lows : 
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Chapter  4 


STRAIN-DISPLACEMENT  RELATIONS 


The  general  nonlinear  strain-displacement  relations  for  large 
rotation  but  small  strain  were  derived  by  Novozhilov  [18]  and  later 
corrected  by  Tsao  [19].  For  shells  of  revolution  with  axisymmetric 
loading,  the  strain  displacement  relations  can  be  written  as 


X 

c 


du 
c 

dS 


dw 


c 


cos  $ + w 

u 


c 


c c 

dS  " R 


sin  <t> ) 


(16) 


The  strains  defined  in  Equations  (16)  are  now  transformed  into 

dx 

the  local  rectilinear  coordinates  as  follows  (recall  that  dS  = cog  g 


and  p = P^ 


du  ~ dw  12 

e,  = cos  p + — cos  P sin  P + 7(Y„) 

1 dx  dx  Z c 


e„  = — (u  sin  V + w cos  ) 
2 r r r 


12 


13 


Chapter  5 


SHELL  DISPLACEMENTS 


The  shell  displacements  are  written  in  local  rectilinear  coordi- 
nates. They  are  represented  by  four  degrees  of  freedom  at  a node:  two 
translations,  one  bending  rotation  and  one  transverse  shear  rotation. 


u = u cos  6 + w sin  p 
c r r 

w = -u  sin  £3  + w cos  (3 
c r r 

dw  / dw  \ , 

_E.  = (_c)dx  (1{ 

dS  \ dx  / dS 

dx  / dur  . R , dWr  R d£  . . dg\ 

= dS  V “dT  Sin  P IT  C°S  P " Ur  COS  P dx  " Wr  Sln  P tej 

(du  dw  u w \ 

- sin  P + cos  P + ~ ^ tan  P ) 

As  was  done  by  Khojasteh-Bakht  [10],  the  displacement  field  is 
assumed  to  be  represented  by 


u = a.  + Q„x 
r I 2 


2 3 

w = Ql  + a,  x + Qcx  + ax 
r 3 4 5 b 


dw  u 
c _c 

r dS  ' r 
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X = cos  p 

r 


/ 

l dx 


dw 

sin  P + — r-^-  cos  p + 
dx 


u w \ 

— + — tan  p) 
r r / 


du 


a cos  6 w sin  8 
r r 


dw 


= - sin  P cos  (3  + cos  6 + 

ax  *'r 


dx 


2 2 

X * sin  P cos  p + (a^  + 2 Q^x  + 3 Q^x  ) cos  P 


rr  = q7  + V 


In  matrix  notation  this  can  be  written  as: 


(19) 


u 

'i 

X 

0 

0 

0 

0 

0 

o' 

\ 

r 

w 

r 

0 

0 

1 

X 

2 

X 

3 

X 

0 

0 

k 

X 

r 

0 

-sc 

0 

2 

c 

2 

2xc 

3x2c2 

0 

0 

y 

r 

0 

0 

0 

0 

0 

0 

1 

X 

J 

(20) 


where 

s = sin  p,  c = cos  P 
This  can  be  written  symbolically  as 

If  ) = [4>]  {cj 


(21) 


where  id)  is  the  generalized  displacements  vector  for  the  curved  shell 
element.  This  displacement  function  allows  rigid  body  motion  without 


inducing  strains. 
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The  shell  displacements  shown  in  Equation  (20)  represent  4 degrees 
of  freedom  at  a node,  twc  translation,  two  rotation.  The  8 degrees  of 
freedom  connected  with  the  nodes  of  the  element  are  written  as  the 
displacement  vector 


l5rJ  - [url,  wrl,  rl,  7rl,  ur2.  wr2,  r2,  7r2l 


(22) 


where  subscripts  1 and  2 refer  to  the  I and  J nodes  respectively. 

The  generalized  displacements  {cj  are  related  to  the  nodal  point 
displacement  vector  by 


(a)  = [A  ] (5e) 


(23) 


ia)  is  evaluated  as  follows: 


1)  At  x = 0 


a,  = u . 
1 rl 


Q~  = w . 
3 rl 


a7  = ’'rl 


-C^  sin  P^  cos  P^  + cos  p^  = 


2)  At  x = l 


U1  + V = Ur2 


2 3 

OL  + a..  I + a.  A + Ol,1  = w 0 

3 4 j o 


+ aa£  - 7 


°7  + a8 


r2 


2 2 2 2 

- 0^  sin  P2  cos  P2  + cos  P2  + 2a,.  i!  cos  P2  + 3a^£  cos  P2  = Xr2< 


Solving  this  system  of  equations  gives 
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where 


19 


20 


where 


|A1  » 


Substituting  (25)  into  (23)  gives 


{cJ 

n 

> 

(6  }e 

r 

“ [Ar] 

[R]  (&}e 

*=  [A]  {&}6 

[A] 

- [Ar] 

[R] 

sin  v 

-cos  i 

0 

0 

0 

0 

0 

0 

-sin  , 

COS  V 

0 

0 

sin  ♦ 

-cos  V 

0 

0 

• 

p 

l 

l 

COS  v 

sin  u 

0 

0 

0 

0 

0 

0 

-a^  sin  , 

cos  * 

bl 

0 

a1  sin  i* 

-a^  cos  * 

0 

0 

a9  sin  v 

-a.,  cos  '4- 

-a2  sin  i|> 

a^  cos  v 

-3  cos  , 

-3  sin  4 

-2b2 

0 

+ 3 cos  v 

+ 3 sin 

K 

0 

.2 

l 

.2 

22 

£2 

b4 

-a^  sin  » 

a ^ cos  ♦ 

a3  sin  v 

-a^  cos  v 

+ 2 cos  , 

+ 2 sin  . 

0 

- 2 cos  ■» 

- 2 sin  , 

K 

Q 

,3 

.3 

b3 

,3 

23 

b5 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

t 

0 

0 

0 

1 

t 

(27) 


(28) 


Chapter  6 


STRESS-STRAIN  RELATIONS 


For  an  axisymmetr ic  shell  of 
loadings,  the  stress  resultants  and 


N1 

'Ell 

E12 

E13 

N2 

E21 

E22 

E23 

M1 

' = 

E31 

E32 

E33 

M2 

E41 

E42 

E43 

<1 

0 

0 

0 

revolution  subjected  to  axisymmetric 
couples  can  be  expressed  as 


E14  0 ' 

el 

E24  ° 

e2 

E34  0 

*1 

E44  0 

< 

2 

0 V 

,'7l 

The  quantities  are  related  to  the  principal  curvilinear  coordi- 
nate system  with  1 as  meridional  direction  and  2 as  circumferential 
direction.  Symbolically  this  can  be  written  as 

(s)  - [E]  (e)  (30) 

where  [E]  is  the  elasticity  matrix.  The  detail  derivation  of  [E]  is 
given  in  Appendix  A. 

Substituting  (19)  into  (17)  gives 

(e)  - [O']  {cJ 

- [O' ] [A]  {6>e  - [B]  (&)e  (31) 


21 


22 


[B]  = [*'][A] 

Differentiating  the  displacement  functions  with  respect  to  x gives 


dx  “ a2 


if  ° °4  + 2Q5  X + 3U6  ^ 


— - * a 

dx  8 


T " 2a5  + 6U6  X 


Substituting  these  relations  into  equation  (17)  gives 


el  = a2  cos  P + (a4  + 2a5  * + 3o ^ x ) cos  p sin  p + j(Xc>2 

e2  = r [^Ql  + Q2  x)  sin  ^ + (a3  + Q4  x + a5  x2  + Q6  x3)  cos  \|>] 

Kl  = -C^3P  (2Q5  + 6Q6  x)  - ^°-s  P sin  P (a4  + 2a5  x + 3a6  x2) 


2 2 

. sin  P - cos  6 . . 

+ ^ H (a2 ) - Qg  cos  p 


*2  = - 


^ [cos2  P («4  + 2a5  x + 3afi  x2)  - sin  P cos  p a2J 


cos  <t>  . . 

~r~  <a7 + n8  x) 


t1  * -(cij  + Qg  x) 


r 


r 
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From  these  equations,  the  matrix  [O' I may  be  obtained  and  may  be 
split  into  two  matrices  [$,£]  and  [4>,n]  containing  linear  and  nonlinear 
terms . 


U)  - [*•]  {a} 

- + l*’"])  (a) 


where  [ <t  ’ ^ ] is  given  by 


■ in  t*  3 k itn  r to* 


co*  , m CO*  | 


-2  co»  . 


co»  : sin  ? -ft  : 


-3  x cos  : co» 


0 0 


and  the  non-zero  first  row  of  [O'11]  is 

2 


(32) 


[®'n]  = 


-Sinp_cosJx  0 copy  xcos2pX  |x2cos2pX  0 0 


(32a) 


The  subscript  1 is  used  to  denote  the  first  row  of  the  matrix. 
From  Equations  (13)  and  (21) 


If  } 

r 


if} 

{f} 


M (a)  - [o>]  [a]  (&}e 
[qrf 1 tfr>  = tqrlT  (tj 

(qrlT  [♦!  [A]  Is}6 
[N]  (&)e 


(33) 


J 
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where 

INI  - [qrlT  [*]  [A]  . (34) 

The  element  stiffness  matrix  and  equivalent  nodal  force  may  be 
obtained  from  the  following  formulas: 


[k6] 

lFe) 

P 


J7  [B]T  [E]  [B]  dA 

A 

e 

//  [N]T  (p)  dA 
A 


(35) 

(36) 


where  (pj  is  the  surface  traction  vector.  The  derivation  of  Equations 
(35)  and  (36)  is  given  in  Appendix  B. 


1 


x 


w 


T 


0 

0 

0 

0 

0 

0 


0 

0 

1 

x 

2 

x 

3 

x 

0 

0 


0 


-sc 


0 

2 

c 


2 

3 

0 

0 


2 

x c 


2 2 


X c 


0 

0 

0 

0 

0 

0 

1 

x 


(37) 


sin  ^ 

-cos  + 

0 

0 

cos 

sin  41 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

(38) 


p = p cos  8 - p sin  8 

*x  rt  rn 

p = p.  sin  8 + P cos  P 

ry  t rn 


prJ  = I qr 1 p) 
p)  * hrlT  Pr) 

sin  i>  cos  t 0 

-cos  it  sin  . 0 

0 0 1 

0 0 0 

Px  sin  , + cos  . 

-Px  cos  v + Py  sin  w 

0 

0 


- 

0 

0 


(41) 

(42) 

(43) 

(44) 


L 
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2 2 
p sin  , + p cos  , sin  v + p cos  i - p sin  \j  cos  v 
x ry  rx  ry 

. 2 . 2 
Px  x sin  , -f  p x cos  i,  sin  t + p x cos  \,  - py  x cos  v sin 

2 2 
p cos  , sin  » + p cos  , - p cos  v sin  v + p sin  ; 


This  value  may  now  be  substituted  into  equations  (34)  and 

(36)  to  obtain  \.Fe}. 

P 


|q  1 ip)  =1 


2 , 2 
p x cos  » sin  . + p x cos  y,  - p x sin  , cos  i + p sin 
x y Kx  Ky 

2 2 2 2 2 2 

p x sin  , cos  * + p x cos  y.  - p x sin  . cos  , + p x sin 

x y rx  v 

3 3 2 3 3 2 

p x sin  . cos  , + p x cos  . - p x s m , cos  , + p x sin 

x rv  r x rv 


Chapter  7 


RESULTS  AND  CONCLUSIONS 

NUMERICAL  EXAMPLES 

To  demonstrate  the  numerical  accuracy  of  the  method,  some  selected 
problems  were  solved  and  compared  to  known  results. 

First,  a circular,  monolithic,  thin  plate  with  clamped  edges  and 
subjected  to  a uniformly  distributed  load  was  considered.  The  plate 
had  a radius  of  100  inches  and  a thickness  of  1,0  inch.  Young's 
modulus  was  1 X 10^  psi  and  Poisson's  ratio  was  0.3.  The  plate  was 
divided  into  five  elements. 

Using  only  five  elements , the  results  from  this  program  agree  to 
within  77o  of  the  exact  results  using  large  deflect  ion  theory  shown  in 
[20] . For  the  particular  loading  case,  elementary  theory  differs  from 
the  exact  solution  by  over  70%.  The  results  are  shown  in  Figure  2. 

A comparison  with  Klein's  [12]  solution  using  linear  analysis 
is  shown  in  Figure  3.  This  was  the  analysis  of  a circular,  flat  plate 
under  axisymmetric  pressure  loading.  It  can  be  seen  that  as  the  number 
of  elements  increased  in  the  linear  solution,  it  approached  the  nonlinear 
solution  using  only  five  elements. 

A hemispherical  shell  is  shown  in  Figure  4.  A comparison  was 
made  with  the  theoretical  values  presented  in  [2],  The  values  calcu- 
lated using  the  0RTH02  program  agreed  with  the  theoretical  values 
within  the  accuracy  with  which  the  curves  could  be  read. 
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HORIZONTAL  DISPLACEMENT  x 10  3 in 
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NOOAL  SPACING 


ANGLE  BETWEEN  NORMAL  TO  SURFACE  AND 
AXIS  OF  ROTATION  (d«g) 


Figure  4.  A Hemispherical  Shell  Solution  by  Finite  Elements 
(Grafton  and  Strome,  J.A.I.A.A.,  1963). 


The  radial  deflection  and  meridional  moment  for  a cylindrical 
shell  subjected  to  a unit  edge  load  is  shown  in  Figures  5 and  6.  It 
can  be  seen  that  the  present  0RTH02  program  agrees  very  closely  with  the 
exact  solution.  These  solutions  are  for  a linear  analysis. 


1 lb /in. 


E - 10  x 10°  psi 
v -0.30 


Figure  5.  Cylindrical  Shell  with  Unit  Edge  Load. 
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x (in.) 

Figure  6.  Meridional  Moment  for  the 
Cylindrical  Shell. 


A spherical  shell  under  uniform  normal  pressure  was  analyzed. 

The  shell  and  its  properties  are  shown  in  Figure  7.  A comparison  was 
made  with  the  exact  solution  given  by  [18].  It  is  readily  seen  that  the 
present  solution  agrees  with  the  exact  solution  very  well. 

The  numerical  influence  of  the  shear  deformation  becomes  much 
clearer  when  a circular  sandwich  plate  with  clamped  edges  subjected  to 
a distributed  load  of  14  psi  is  considered.  The  plate  has  a radius  of 
10  in.,  the  thickness  of  core  layer  is  0.75  in.,  and  the  thickness  of 
upper  and  lower  facings  is  0.028  in.  and  0.022  in.,  respectively.  Young's 
modulus  of  facings  is  10^  psi,  Poisson's  ration  is  0.3,  and  the  shear 
modulus  of  core  is  30,000  psi.  The  plate  is  also  divided  into  5,  10, 
and  20  elements.  The  results  are  given  in  Figure  8.  The  maximum  deflec- 
tion of  the  plate  is  shown  to  converge  to  the  theoretical  value  of 
0.0415  in.  as  reported  by  Plantema  [19], 


EXACT 


Figure  8.  Convergence  of  Center  Deflection 
of  Circular  Sandwich  Plate. 

Nickell  [14]  obtained  a solution  for  a cylinder  loaded  with  a 
radial  load  on  one  end  and  a moment  (see  Figure  9),  The  results  are 
compared  with  the  present  solutions  in  Figures  10  through  12.  The 
results  agree  with  Nickell's  within  the  limits  of  accuracy  with  which 
the  curves  can  be  read. 


R > 18  5 in 
t ■ 3.0  in. 

L ■ 36.0  in. 

E - 3*  106  in. 

i’  * 0 

P - 1500  lb /in. 

M - 1000  lb— in. /in. 


Figure  9.  Locally  Loaded  Cylinder. 


Figure  12.  Shear  Versus  Length. 


Sharifi  [13]  analyzed  a clamped  circular  sandwich  plate  under 
a uniform  lateral  pressure.  A comparison  with  his  solution  and  the 
linear  solution  is  shown  in  Figure  13.  He  used  an  incremental  formu- 
lation for  a nonlinear  finite  element  analysis  of  sandwich  structures. 
The  nonlinearities  considered  were  due  to  large  displacements,  as  a 
result  of  finite  rotations,  and  plastic  deformations  of  the  facings. 
The  0RTH02  program  solution  agreed  very  closely  with  these  results 
as  evidenced  in  Figure  13.  The  nonlinear  solution  differs  signifi- 


cantly from  the 


linear . 
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A finite  element  method  that  used  the  displacement  model  was 
selected  to  analyze  the  system.  The  meridional  curve  of  the  shell  was 
represented  by  a series  of  curved  elements  having  local  coordinates. 

An  element  was  developed  that  matches  slopes  and  curvatures  as  well  as 
displacements  at  its  nodal  circle. 

The  geometric  approximation  of  curved  shells  usually  associated 
with  the  finite  element  method  is  minimized  by  using  the  curved  element. 
The  use  of  this  curved  element  significantly  reduces  the  meridional 
bending  moment  usually  present  at  the  nodal  circles  when  a curved  struc- 
ture is  approximated  by  a straight  line  (conical)  segment.  A smaller 
number  of  elements  can  be  used  in  comparison  to  that  of  a conical 
e lement . 

A computer  program  was  developed  to  solve  the  derived  equations. 
The  program  was  shown  to  be  a versatile  and  flexible  method  of  imple- 
menting the  basic  theory.  Several  problems  were  solved  and  the  results 
compared  with  both  linear  and  nonlinear  solutions  from  the  literature. 

In  most  cases,  the  results  from  this  program  agreed  with  the  linear 
solutions  within  the  limits  of  accuracy  with  which  the  curves  could  be 
read.  Agreement  with  nonlinear  solutions  was  good  and  could  usually  be 
further  improved  by  taking  more  elements.  The  shell  thickness  and 
pressure  may  vary  linearly  along  the  meridian.  The  convergence  and 
accuracy  of  the  method  were  found  to  be  entirely  satisfactory  as  evi- 
denced by  the  numerical  examples. 

The  Gaussian  Quadrature  Integration  method  was  used  in  the 
derivation  of  the  stiffness  matrix.  Several  tests  were  made  to  determine 
the  most  efficient  method.  As  many  as  eight  points  were  used.  After 
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examining  the  different  schemes,  it  was  decided  that  the  two  point 
Gaussian  Integration  scheme  gave  the  best  results. 

The  accuracy  obtained  by  this  method  depends  directly  on  the 
extent  to  which  the  assumed  displacement  patterns  are  able  to  reproduce 
the  deformation  actually  developed  within  the  element.  Since  the  chosen 
displacement  patterns  satisfy  the  requirements  of  completeness  and 
conformity  (continuity  of  displacement  at  element  boundary)  as  the  size 
of  the  element  decreases  indefinitely,  the  solution  obtained  converges 
to  the  exact  solution. 

The  finite  element  method  is  obviously  a powerful  tool  in  the 
analysis  of  orthotropic,  composite  structures.  There  still  remains 
much  work  to  do  in  this  area.  A logical  extension  of  this  work  is  to 
include  stability  criteria.  Other  items  which  should  be  considered  in 
the  future  are:  crossover  effects,  cracking  or  "crazing"  of  the  matrix 
material,  an  appropriate  failure  criterion,  and  material  properties 
which  are  different  in  tension  and  compression. 
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Appendix  A 

ELASTICITY  MATRIX 

Individual  curved  finite  elements  can,  in  general,  be  composed 
of  a number  of  anisotropic  layers  of  varying  thickness  along  the 
meridional  coordinate.  For  a single  lamina,  considering  shear  deforma- 
tions, the  constitutive  relation  is  given  as 
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where  the  transverse  normal  stress  has  been  omitted  and  the  laminae 

L} 

are  orthotropic  with  respect  to  the  principal  elastic  axes  L-T.  These 
axes  need  not  coincide  with  the  axes  of  the  curvilinear  coordinate  system 
1-2,  (Figure  14),  (1  is  the  meridional  direction)  and 
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To  develop  a theory  for  structural  laminates  with  individual 
layers  having  their  elastic  axes  oriented  at  various  angles  relative  to 
the  coordinate  axes,  the  stress-strain  Equations  (A-3)  must  be  rotated 
through  the  positive  angle  0 so  that  the  transformed  stress-strain 
equations  are 


and 
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Substituting  the  midsurface  strain  and  curvatures  into  Equations 
(A-4)  the  following  expression  is  obtained: 


and 


By  integrating  over  the  total  thickness  of  the  laminate,  the 
generalized  stress  resultants  in  terms  of  midsurface  strain  and  curva- 
ture are  given  as 
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where 


[C]  = Z [Q  (1°]  (h  - h ) 
k=l 


IB-1  ' I S [Q(k)l  (h2k  - K2.,) 


k=  1 


[D] 


- i “ (Q(k)i  <hk  - hk-i} 


k=l 


[S]  = Z tQ(k)l  (h  - hkl> 
k=l 


(A-8) 


in  which  h^  and  h^  ^ = the  distances,  respectively,  from  the  midsurface 
to  the  inner  and  outer  surfaces  of  the  k-th  layer. 

For  an  axisymmetric  shell  of  revolution  subjected  to  axisymmetric 
loadings,  N12  - M12  = Q2  * el2  = 'c12  * y 2^  = °* 

Hence, 
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(s)  - [E]  {e}  . (A-10) 


or  symbolically 
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Appendix  B 

ELEMENT  STIFFNESS  MATRIX 

The  element  stiffness  matrix  is  found  by  writing  the  total 
potential  energy  of  the  axisymmetric  shell  of  revolution  and  minimizing 
it  for  the  imposed  constraints  and  loading  conditions. 

The  potential  energy  for  a linear  elastic  shell  of  revolution 
in  the  absence  of  thermal  and  body  forces  can  be  formulated  as  follows: 

it  = ///  \ (e)T  (a)  dV  - //  (f)T  (p)  dA  (B-l) 

V A1 

where  the  vectors  (e),  (cr),  If),  and  (p)  represent  the  strain,  stress, 
displacement,  and  equivalent  surface  traction  vectors,  respectively. 
Introducing  the  stress  resultant  vector 

is)  - t (a)  (B-2 ) 

where  t is  the  thickness  of  the  shell,  Equation  (B-l)  may  be  written  as 

it  = ///  \ WT  (s)  f - a {f)T  (p)  dA  . (B-3) 

V A1 

The  first  integral  is  evaluated  over  the  entire  volume  V of  the 
shell  and  the  second  over  the  portion  A^  of  the  midsurface  of  the  shell, 
where  the  equivalent  surface  tractions  are  prescribed.  Since  the  state 
of  displacement  throughout  the  shell  is  defined  element  by  element,  the 
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total  potential  energy  may  be  considered  as  the  sum  of  the  potential 
energies  of  all  individual  elements,  i.e., 


t 


The  potential  energy  contribution  of  element  "e"  will  now  be 
considered.  The  state  of  displacement  defined  for  the  element  in  local 
rectilinear  coordinates  x-y  can  be  expressed  in  matrix  form  in  Equation 
(21)  as 


{frl  = [0]  (ct)  = [<t>][Ar]  lb®)  . (B-4) 

Transformation  of  {ff}  into  the  global  coordinate  system  may  be 
obtained  from  Equation  (13) 

If)  = [qr]T  {fr)  = [N]  (5r}e  (B-5 ) 

where 

[N]  = [qr]T  [0]  [Ar]  (B-6) 

and  the  colume  vector  represents  the  eight  discrete  parameters 

(nodal  point  displacements)  of  the  element  as  given  in  Equation  (25b). 
The  matrix  [N]  is  a function  of  spatial  coordinates  and  describes  the 
defined  displacement  pattern. 

Substituting  Equation  (27)  into  Equation  (31)  the  following 


strain-displacement  relations  are  obtained: 

(e)  - [B]  (&}e 


(B-7) 
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where 

[B]  = [♦']  [A]  . (B-8) 

Equation  (B-8)  is  a matrix  relating  the  nodal  point  displacement  vector 
to  the  strain  vector.  The  elastic  stress-strain  relations  can  be 
expressed  as 

(s)  = [E]  (e)  (B-9 ) 

where  [E]  is  a function  of  the  elastic  properties  of  the  element.  Each 
element  can  be  assigned  different  elastic  properties.  If  the  relations 
in  Equations  (B-9),  (B-5),  and  (B-7)  are  substituted  into  (B-3),  the 
potential  energy  contribution  for  the  element  becomes 

*e  “ fff  \ [B]T  [EJ  [B]  (6e)  Y - //  tNlT  ^ ^ 

Ve  A1  (B- 10) 

e 


where  V is  the  volume  of  the  element  and  A,  is  that  part  of  the  mid- 
e 1 

e 

surface  area  of  the  element  which  coincides  with  the  midsurface  area  A^ 
of  the  shell  over  which  the  equivalent  surface  tractions  are  prescribed. 

Since  the  discrete  parameters  (&ej  are  not  a function  of  spatial 
coordinates,  the  potential  energy  of  the  element  may  be  written  as 


e 


n 


fff  J [b]t 

V 


[E][B] 


lb6}  - (&e}T  //  [N]T  (p)  dA  . 

le  (B-ll) 


Since  the  assumed  displacement  patterns  for  each  element  satisfy 
various  requirements  such  as  completeness  and  conformity,  the  best 
values  that  can  be  obtained  for  the  total  nodal  point  displacements  of 
the  finite  element  representation  of  shells  of  revolution  are  those 


that  minimize  the  total  potential  energy  of  the  shell  under  the  con- 
straints imposed;  i.e.,  the  best  value  of  {b}  are  those  that  satisfy 
the  system  of  linear  equations 
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-2-  = 0 


(B-12) 


where  ib)  is  the  total  nodal  displacement  vector  of  the  system. 

in  forming  the  system  of  Equations  (B-12),  it  is  convenient  to 
have  an  expression  for  the  spatial  derivatives  of  the  potential  energy 
of  each  element  "e"  with  respect  to  its  own  nodal  point  displacement 
vector  lbej , i.e . , 


x e e^e  e > e e e e 
on  on  on  on  on  on  on  on 


sibe)  Lcui  L'xi  'xj 


(B-13) 


By  use  of  Equation  (B-10),  this  expression  can  be  obtained  as 


5ne 

d(be) 


///  [B]1  [ E ] [ B ] 

V 


{&e}  -I" / J 

J 


//  [N]1  Ip/  dA  . (B-14) 


The  terms  in  the  first  and  second  brackets  are  normally  defined 

0 

as  the  element  stiffness  matrix  [K  ] and  the  element  generalized  nodal 
point  force  (f6),  respectively.  Hence, 


[Kel  = ///  t B ] T [ E]  [ B]  ^ 

V C 

e 


(B- 15  ) 


If6}  - //  [n]t  (p)  dA 
A, 


(B-16) 
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By  properly  combining  the  submatrices  in  Equation  (B-14)  obtained 
for  each  element,  the  total  matrix  equation  representing  Equation  (B-12) 
can  be  constructed  as 

[K]  l&}  = {fI  (B-17) 

and  then  solved  for  the  nodal  point  displacements.  Once  the  nodal 
point  displacements  are  obtained,  the  corresponding  stress  resultants, 
stresses,  and  strains  for  the  defined  displacement  patterns  can  be 
calculated  from  Equations  (B-7)  and  (B-9). 

If  Equation  (B-8)  is  substituted  into  Equation  (B-15)  and  the 
volume  increment  for  a shell  of  revolution  is  taken  as 

dV  = 2nt  d*  , (B-18) 

cos  fJ 

then  the  element  stiffness  matrix  for  the  axisymmetric  shell  element 
takes  the  form 


[Ke]  = 2k  f*  [B]T  [E]  [B]  dx 

q cos  p 

= 2k  [A]T  [G] [A]  (B-19) 

where 

[G]  = / [O’  ]T  [E][*’]  -^r  dx  . (B-20) 

o cos  P 

The  integration  is  over  the  chord  length  of  the  meridian  cross 
section  of  element. 

It  is  assumed  that  the  equivalent  surface  traction  over  the  mid- 
surface area  where  tractions  are  prescribed  varies  linearly  between 
the  two  nodal  circles  I and  J.  That  is, 


(p  }T  = [0  (P  + p'  X)  0 0] 


(B-2 1 ) 


where  (p  ) is  the  surface  traction  vector  expressed  in  local  curvilinear 
c 

coordinates.  Transforming  into  global  coordinates  the  following  is 
obtained : 


(p)  = [qr]T  tqclT  (p<J 


(B-22) 


Substituting  Equations  (B-6)  and  (B-22)  into  Equation  (B-16)  the 
generalized  element  nodal  force  vector  becomes 


where 


(Fe)  = 2 it  / [ArlT  [^1T  [qrU^r]T  hclT  CRJS*  ft  dx  (B-23 


hV  - 2»iAriT/  mT  lic)T  tPc-l  d* 


-sin  p 


-x  sin  p 


mT  [qjT  iO 


-x  sin  P 
cos  P 


x cos  P 


x cos  P 


-x  sin  P 
x cos  P 


x cos  P 


x cos  P 


x cos  P 
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Appendix  C 


EQUILIBRIUM  EQUATIONS  BY  VIRTUAL  WORK 

The  equilibrium  equations  governing  the  shell  behavior  can  be 
derived  by  using  the  principle  of  virtual  work. 

From  Equations  (21),  (23),  and  (25a) 

if)  = [•]  (cj 

= [0]  [A]  (u) 

= [» 1 [A] [R]  (ug)  (C-l) 

“ [ $ ] [A  ] (u  ) 

8 

= [N]  lu  } 
g 

where  {N)  = [i>][A}  and  [A]  = fA][R]  . 

Let  there  be  an  arbitrary  and  non-zero  virtual  nodal  displacement 
&{ugj  about  the  deformed  position  which  results  in  a virtual  displace- 
ment &{ f j and  virtual  strain  bit).  The  & prefix  denotes  a virtual 
change  in  the  quantity  concerned. 

By  means  of  the  principle  of  virtual  work,  the  following  expression 
can  be  written: 

6(u  )T  (q)  » / bit}1  (t)  dA  - / If  ) dA  (C-2) 

g i m A S 

m m 

where  A is  the  reference  surface  area  of  the  shell,  {q}  is  the  applied 
m 

nodal  force  vector  and  (f  ) is  the  surface  traction  vector.  The  stress 

s 
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vector  (t)  is  given  as 


(t)  - [C ] (e) 


[C][*'][A]{u  } 


Substituting  Equation  (C-l)  into  (C-2)  the  following  is  obtained: 


&lu  )T  1q)  + &lu  )T  / [N]T  (f  ) dA  « / 6(e}T  {T} 


&(u  )T  (p)  = / 5(e}T  { t} 

o J 


(C-3) 


where  { P j is  the  equivalent  nodal  forces  of  the  element  defined  by  the 
principle  of  virtual  work. 

Substituting  Equations  (23),  (25a),  and  (31)  into  (C-3)  gives: 


&lu  )T  (p)  - / &(u  )T  [A  ] C [0']T  It) 


(C-4) 


This  results  in  a nonlinear  matrix  equation  for  the  equivalent 
nodal  forces  (p),  i.e.. 


(p)  - / (A]T  [O' ]T  [C]  [O' ] [A]  (u  } dA 
A 8 m 


(C-5  ) 


Equation  (C-5)  is  now  linearized  by  writing  it  in  the  form  of  an 
implicit  differential. 
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A {p>  ■ / A [A]T  [f ' ]T  tCHO'HAl  {u  ) dA 

A 8 m 

m 

+ / [A]T  a(i«'1T  ])  [A]  in  } dA  (C-6 ) 

Am  V y g ra 

+ / (A]T  [♦* ] |C] [•• ] [A]  A (u  } dA 

A 8 m 

m 


It  is  assumed  that  a change  in  the  transformation  matrix  during 
an  increment  of  load  may  be  neglected.  This  permits  neglecting  the 
first  term  on  the  right  hand  side  of  equation  (C-6).  The  second  term 
results  in  the  well  known  initial  stress  matrix  while  the  third  term 
accounts  for  the  effect  of  the  increment  of  strain  and  may  be  split 
into  two  terms  separating  the  linear  and  nonlinear  displacement  terms. 
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(C-7) 
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(2 ) 

The  last  three  terms  have  been  collected  into  [ K ] and  will  be 
called  the  initial  displacement  matrix.  If  the  area  increment  of  the 
shell  is  taken  as 


then  the  element 
takes  the  form 


dA 


2rt 


3JjlL 


dx 


cos  P 

stiffness  matrix  for  the  axisymmetric 


shell  element 


[Ke]  - 2ji  / [B]T  [C]  [B]  R (x)  dx 
o 

- 2n  [A]T  [G] [A] 

where 

[B]  - 1 [A] 

and 

[G]  = f [*’ ]T  [C] [0 • ] R(x)  dx  . 
o 

The  integration  is  over  the  chord  length  of  the  meridian  cross 
section  of  element. 


INITIAL  STRESS  MATRIX 


The  second  term  of  Equation  (C-6)  can  be  written  as: 
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2/ 
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2/ 
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[A]T  a([4>']T  t C]  [ 0> ' ])  [A]  (u  } dA  « 
V / g m 

[A]T  (a[*']T)  [c][0'  ] [A]  iu  }dA  « 
[A]T  ^AI  «>  ’ ] T ) if)  dAm  . 


(C-8) 
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Equation  (C-8)  can  be  broken  down  into  a summation  of  the  five 
stress  resultant  components,  i.e.. 


2/  [A]1  A[*'  ] 1 (t)  dA  - £ / t.[A]T  2A  {$'.} 

» m . , * 1 l 

A l-l  A 

m m 


where  i is  the  index  of  the  stress  resultant  components  and  is 

the  transpose  of  the  corresponding  i1"*1  row  of  the  matrix  [ 4> ' ] . 

Since  the  derivation  is  based  on  the  current  deformed  position 
of  the  shell  element,  the  increment  aI'I’!^'}  can  be  written  as  follows; 


2a(  $ ' T)  = 2a{  <t>  J nTj  = <»• ! nT> 

it  1 


and 
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Appendix  D 


THEORETICAL  BACKGROUND 

Some  of  the  basic  concepts  used  in  the  derivation  of  the  equations 
presented  in  the  main  text  are  presented  in  this  section. 

GEOMETRY  OF  SHELLS 

The  geometry  of  a shell  is  entirely  defined  once  the  midsurface 
and  the  thickness  at  each  point  are  specified . Hence,  to  describe  the 
shell  space,  the  middle  surface  or  reference  surface  of  the  shell  must 
be  specified.  Let  and  be  the  curvilinear  coordinates  for  the 
mid-surface  and  let  them  coincide  with  lines  of  principal  curvature  of 
the  surface,  and  let  £ be  a coordinate  normal  to  the  midsurface  as  shown 
in  Figure  15. 


Figure  15.  Typical  Shell  Element. 
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The  position  of  any  point  in  the  midsurface  can  be  defined  by 
the  curvilinear  coordinates  0^  and  . The  location  of  any  point  in 
the  shell  can  be  related  by  the  three  parameters  , and  With 

this  curvilinear  coordinate  system,  a line  element  in  the  shell  space 
surrounding  the  midsurface  can  be  expressed  in  terms  of  the  differen- 
tials of  the  orthogonal  curvilinear  coordinates  as  follows: 


where  A^  and  are  the  midsurface  metrics  and  and  are  the 
principal  radii  of  curvature  of  the  surface. 


STRAIN- DISPLACEMENT  RELATIONSHIPS 

The  general  nonlinear  strain-displacement  relations  for  large 
rotation  but  small  strain  were  derived  by  Novozhilov  [18]  and  later 
corrected  by  Tsao  [19].  Suppressing  the  nonlinear  terms,  the  following 
strain-displacement  relations  for  linear  theory  of  shells  is  obtained: 
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where  the  functions  U,  V,  and  W represent  the  displacement  components 
caused  by  straining  of  a material  point  originally  at  point  (G , 
in  the  shell  in  the  G^ , Q^,  and  £ direction,  respectively 

To  incorporate  the  transverse  shear  deformation,  the  classical 
Kirchhof f-Love  assumption  must  be  abandoned.  The  material  lines  origin- 
ally straight  and  normal  to  the  midsurface  of  the  shell  remain  straight 
but  are  no  longer  normal  to  the  deformed  midsurface  (Figure  16).  This 
implies  that  the  transverse  shear  deformation  is  independent  of  the 
coordinate  £.  Hence,  the  shear  rotation  can  be  represented  by  some 
average  value  of  the  shear  strain  at  midsurface.  The  displacement 
components  of  a point  in  the  shell  can  be  expressed,  as  a first  approxi- 
mation, by  relationships  of  the  form 


U(«1,a2,0  = u(alta2)  + t; 

v(alfa2,0  = V(alta2)  + £ 
w(a1,a2,0  - w(a1,a2) 


(D-3) 


where  u,  v,  and  w are  displacements  of  the  point  at  midsurface,  and  (3^ 
and  P2  are  rotations  that  represent  changes  of  slope  of  the  normal  to 
the  midsurface.  It  should  be  noted  that  terms  u,  v,  w,  and  P2  are 
functions  of  and  only. 


UNDEFORMED  STATE 


DEFORMED  STATE \ \ 


Figure  16.  Transverse  Shear  Deformation. 


Substituting  Equations  (D-3)  into  the  strain-displacement  rela- 
tions Equation  (D-2)  and  suppressing  the  terms  yield 
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are  the  extensional  strains  at  the  midsurface  of  the  shell, 
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are  the  changes  in  curvature  of  the  midsurface  in  the  directions  of 
and  u , respectively,  and 
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represent  the  in-surface  shear  strain  and  twist  of  the  midsurface, 
respect ive ly . 


STRESS-STRAIN  RELATIONS 


Assuming  that  the  in-surface  stresses  can  be  represented  by  a 
state  of  generalized  plane  stress,  the  stress-strain  relations  for  the 
space  and  for  the  orthotropic  material  can  be  written  as 
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T 12  = G12  712 
Tl£  = Gl^  71£ 
t2^  = G2 £ 71£ 


(D-8  ) 


where  E^,  E^ , G^.)t  G^,  G , v12  ’ ancl  V21  are  tlle  e^asti-c  constants 
along  the  three  coordinate  directions  [21], 

Since  the  transverse  shear  strain  has  been  assumed  to  be  constant 
across  the  thickness,  the  corresponding  shear  stress  is  likewise  constant 
and  is  directly  proportional  to  the  shear  strain.  However,  from  ele- 
mentary strength  of  materials  it  is  known  that  transverse  shear  stress 
is  not  constant  across  the  thickness  of  a beam  section.  Therefore,  the 
average  shear  strain,  which  may  provide  a good  approximation  to  the 
shear  rotation,  does  not  necessarily  provide  an  adequate  representation 
of  the  shear  stress  distribution.  Hence,  a shear  stress  factor  is  used 
in  conjunction  with  the  transverse  stress-strain  Equation  (D-8)  as 
suggested  by  Naghdi  [22],  that  is 


TiC-|ci^iC  ‘“■2 


(D-9) 


Substituting  Equations  (D-4)  into  Equations  (D-8)  and  (D-9)  and 
integrating  across  the  thickness  of  the  shell,  the  stress  resultants 
and  couples  are  obtained  as  follows: 
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where  t is  the  thickness  of  the  shell  and 
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SHELLS  OF  REVOLUTION 


(D-10) 


(D-ll) 


The  discussion  will  now  be  restricted  to  shells  of  revolution. 
The  midsurface  of  the  shell  is  obtained  by  rotation  of  a plane  curve. 
This  curve  is  called  the  meridian  and  its  plane  is  the  meridian  plane. 
The  intersection  of  the  surface  with  planes  perpendicular  to  the  axis 
of  rotation  are  parallel  circles  and  are  called  parallels.  For  shells 
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of  revolution,  the  lines  of  principal  curvature  are  its  meridians  and 
parallels  [23] . 

A set  of  convected  normal  coordinates  $,  0,  and  £ are  used  to 
describe  the  shells  of  revolution,  where  $ is  the  angle  between  the 
normal  to  the  midsurface  of  the  shell  and  the  axis  of  revolution,  0 is 
the  angle  describing  the  position  of  points  of  the  corresponding  parallel 
as  shoi’n  in  Figure  17.  The  radius  of  curvature  of  the  meridian  is  R^. 

The  second  radius  of  curvature  will  always  be  the  length  of  the 
intercept  of  the  normal  to  the  midsurface  between  the  axis  of  the  shell, 
i.e.,  AP.  This  is  because  the  normal  from  two  adjacent  points  P and  P' 
on  the  same  parallel  will  always  intersect  on  the  axis  of  the  shell. 
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Figure  17.  Shell  Geometry. 

The  arc  length  of  a line  element  in  the  shell  space  is  given  as 
2 2 2 2 2 2 

dS  “ R^  do  + R2  sin  <t>  do  . (D-12) 

Associating  0^  with  <t  and  with  0 and  comparing  Equation  (D-12) 
with  Equation  (D-l),  the  following  is  obtained: 


L 


R = R2  sinO 


(D-13) 


A 


2 


Note  that  the  term  in  Equation  (D-l)  is  small  compared  to  unity  for 

i 

thin  shells. 

From  Figure  18,  by  inspection,  the  following  is  obtained: 


cos<J> 


or 


dR 

d<l> 


cos<t 


(D-14) 


Figure  18.  Shell  Meridian. 

AXISYMMETRIC  LOADINGS 


For  shells  of  revolution  with  axisymmetric  loadings,  all  geometric 
quantities  are  independent  of  0.  Consequently,  all  of  the  shell  vari- 
ables are  independent  of  0 and,  starting  with  the  relationships  between 
the  strains  and  displacements  Equations  (D-4)  through  (D-7),  the 
following  is  obtained: 
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1 du  w 

61  ~ Rx  do  Rl 


e_  *>  — (u  cosO  + w sinO) 
2 R 


1 Rx  dO 


(D-15) 


<2  - y cos<l1 


1 dw  u , Q 

71{;  Rx  dO  R^  ^1 


Setting  y ^ equal  to  some  average  shear  strain  as 


> = "’'I 


(D-16  ) 


and  substituting  Equation  (D-16)  into  the  last  Equation  of  (D-15)  yields 


a - /—  u . \ 

pi  = "yR1  d»  ' rl  7 1 J ’ 


(D-17) 


Furthermore,  let 


X _d i. 

Rx  do  ~ dS 


(D-18) 


where  S is  measured  along  the  meridional  direction  of  the  midsurface. 

Substituting  Equations  (D-18)  and  (D-17)  into  Equation  (D-15)  the 


du  w 
61  = dS  R, 


e0  = — (u  cosO  + w sinO) 
2 R 


following  is  obtained: 
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The  stress  resultants  and  couples  reduce  from  Equation  (D-10)  to  the 
following  set: 


= c, , 

e „ 

+ c,„ 

11 

1 

12 

= c21 

ei 

+ c22 

= D,  , 

#c 

+ D, . 

11 

1 

12 

■ D„ , 

K, 

-f-  D-  n 

21 

1 

22 

-!« 

1 yn 

J 


L 


Appendix  E 


COMPUTER  PROGRAM  0RTH02 

The  foregoing  derivation  was  implemented  in  a finite  element 
computer  program.  For  convenience  in  reference,  this  program 
will  be  referred  to  as  0RTH02 . The  element  stiffness  matrix  was 
formed  by  numerical  integration.  The  nodal  point  coordinate  and  element 
connection  array  are  generated  automatically  by  the  program.  A shell  of 
revolution  is  first  divided  into  as  many  segments  as  necessary.  Because 
each  segment  may  be  considered  as  a separate  unit,  different  material 
properties  as  well  as  thickness  and  pressure  can  be  ascribed  to  different 
segments.  Each  segment  in  turn  may  be  subdivided  into  any  number  of 
shell  elements.  Normal  pressure  and  thickness  of  the  shell  must  be 
axisymmetr ic , but  may  be  varied  linearly  along  the  meridional  direction. 
The  matrix  equations  are  solved  by  the  Cholesky  decomposition  process 
which  stores  only  nonzero  elements  and  therefore  results  in  a significant 
saving  of  computing  time.  The  program  can  be  used  to  solve  problems  of 
thin,  thick,  and  sandwich  shells  of  revolution  as  well  as  multilayered, 
orthotropic  shells  such  as  a fiber  reinforced  composite.  The  program  is 
limited  to  ten  different  materials  and  50  nodes,  but  can  be  increased  by 
increasing  the  dimension  statement  accordingly.  This  program  requires 
about  32K  core  storage  and  three  scratch  files.  A total  of  nine  sets  of 
input  data  is  needed.  The  flow  chart  for  0RTH02  is  shown  in  Figure  19. 

An  iteration  procedure  is  used  in  the  program.  For  the  nonlinear 
effect,  the  load  is  applied  in  increments  and  the  coordinates  are  updated. 
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This  method  was  compared  with  that  of  using  the  relations  in  Equation 
(32a)  and  was  found  to  give  practically  the  same  results.  Much  of  the 
data  is  generated  internally  in  the  program. 


Figure  19 


Flow  Chart  for  0RTH02 . 
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Figure  19 


(Cone luded ) . 
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DATA  INPUT  INSTRUCTIONS 

1.  Problem  card;  Load  Increment  card  (215): 


Col.  1-5 

Number  of  problems  per  run 

(NPROB) 

6-10 

Number  of  load  increments 

(NINCR) 

2. 

Title  card  (9A8): 

Col.  1-72 

Title  to  be  printed  with  output 

(TITLE) 

3. 

Control  card  (515): 

Col.  1-5 

Number  of  boundary  points 

(NBP  ) 

6-10 

Number  of  segments 

(NSEG  ) 

11-15 

Number  of  elements 

(M) 

16-20 

Number  of  nodes 

(NQ1) 

21-25 

Number  of  materials 

(NMAT) 

4. 

Material  cards 

(7F10.0)  one  for  each  material: 

Col.  1-10* 

Young's  modulus  in  meridional  direction 

(El) 

11-20* 

Young's  modulus  in  circumferential 
direction 

(E2) 

21-30 

Poisson's  ratio  in  meridional  direction 

(PR1) 

31-40 

Poisson's  ratio  in  circumferential 
direction 

(PR2 ) 

41-50 

Enter  0.  for  thin  shell 

1.  for  thick  shell  with 

G “ E/2(l+v)  or  shear  modulus 

for  thick  or  sandwich  shell 

(Gl) 

5. 

Boundary  cards 

(515,  5X,  4F10.2)  one  for  each  boundary  point: 

Col.  1-5 

Boundary  node  number 

(I) 

6-10 

r-direction  0 free 

1 fixed 

(ID1) 

11-15 

z-direction  0 free 

1 fixed 

(ID2) 

16-20 

Normal  rotation  0 free 

1 fixed 

(ID3) 

21-25 

Shear  rotation  0 free 

1 fixed 

(ID4) 

31-40 

Prescribed  r displacement 

(UP) 

41-50 

Prescribed  z displacement 

(WP) 

51-60 

Prescribed  rotation 
(angular  displacement) 

(THP  ) 

61-70 

Skewed  boundary  (angle) 

(AL)  j 

L. 
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6.  Segment  cards 

(IX,  212,  15,  5 F10.2,  415)  one  for  each  segment;  1 

Col.  2-3 

0 Conical  segment  (straight  line) 

1 Spherical  segment 

2 Elliptical  segment 

3 Arbitrary  curved  segment 

(ICODE) 

4-5 

Number  of  layers 

(NLAYER) 

6-10 

0 NLAYER  is  the  same  as  the  previous 
segment 

(LAY ID) 

1 New  NLAYER  for  the  segment 

New  layer  data  are  required 

11-20 

R-coordinate  of  the  first  node  of  the 
segment 

(Rl) 

21-30 

Z-coordinate  of  the  first  node  of  the 
segment 

(Zl) 

1 

31-40 

Total  length  of  the  segment  if 

ICODE  - 0 

(Al) 

Total  subtend  angle  of  the  segment  if 
ICODE  «=  1 

The  difference  in  the  R-coordinate  of  the 
first  and  last  node  of  the  segment  if 

ICODE  = 2 

Blank  if  ICODE  = 3 

41-50 

Angle  of  slope  between  the  straight  line 
segment  and  the  r-axis  if  ICODE  = 0 

(A2) 

Radius  of  the  spherical  segment  if 

ICODE  = 1 

Major  radius  of  the  elliptic  segment  if 
ICODE  « 2 

Blank  if  ICODE  = 3 

51-60 

Leave  blank  if  ICODE  = 0 

(A3) 

Phase  angle  <t>  between  the  normal  to  the 
shell  surface  and  the  axis  of  revolution 
if  ICODE  = 1 (to  the  first  node  of  the 
segment ) 

Minor  radius  of  the  elliptic  segment  if 
ICODE  - 2 

Blank  if  ICODE  = 3 

61-65 

First  element  number  of  the  segment 

(Ml) 

66-70 

Last  element  number  of  the  segment 

(M2) 

71-75 

First  node  number  of  the  segment 

(Nl) 

76-80 

Last  node  number  of  the  segment 

(N2) 
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Pressure  loading  and  thickness  cards  (8F10.2)  one  for  each  seg- 
ment except  ICODE  = 3.  (Replace  7a  by  7b  and  7c  when  ICODE  = 3): 


Col.  1-10 

Normal  pressure  at  the  first  node  of 
the  segment 

(PI) 

11-20 

Normal  pressure  at  the  last  node  of 
the  segment 

(P2) 

21-30 

Thickness  of  the  shell  or  thickness  of 
the  core  layer  of  a sandwich  shell  at 
the  first  node  of  the  segment 

(Tl) 

31-40 

Thickness  of  the  shell  or  thickness  of 
the  core  layer  of  a sandwich  shell  at 
the  last  node  of  the  segment. 

(T2) 

Replace  7a  by 

the  following  set  if  ICODE  = 3. 

Element  cards 

(5110)  one  for  each  element: 

Col.  1-10 

Element  number 

(I) 

11-20 

Node  1 1 . 

\ element  connection 

(J) 

21-30 

Node  j) 

(K) 

31-40 

Number  of  layers 

(MLAYER) 

41-50 

Layer  idenfication  code 

(LID) 

Coordinate  cards  (15,  5X,  2F10.2,  2F5.1,  4F10.2)  one 

for  each  node: 

Col.  1-5 

Node  number 

(N3) 

11-20 

R-coordinate  of  the  node 

(R) 

21-30 

Z-coordinate  of  the  node 

(Z) 

31-35 

Angle  between  normal  to  the  shell 
surface  and  the  axis  of  symmetry 

(PHA) 

36-40 

Meridian  curvature  of  the  shell 

(KAPP) 

41-50 

Normal  pressure  at  the  node 

(PP) 

51-60 

Thickness  of  the  shell  or  thickness 
or  the  core  layer  of  a sandwich  shell 

(TT) 

Concentrated  load  cards  (215,  F10.2)**  one  for  each  load 

component  plus 

one  EOD  card: 

Col.  1-5 

Node  number 

(I) 

6-10 

1 for  r-component  of  the  load 

2 for  z-component  of  the  load 

3 for  moment  loading 

(NCI) 

11-20 

Magnitude  of  the  loading.  Positive 
if  the  direction  of  loading  coincides 
with  the  positive  direction  of  the 
coordinates 

(VI) 
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9.  Layer  data  (15,  5X,  3F10.4)  new  set  of  data  is  required  if 
there  is  a change  in  the  number  of  layers: 


1-5 

Material  type  number 

(MT1 ) 

11-20 

Distance  from  reference 
of  layer  at  node  I 

surface 

to  top 

(HU) 

21-30 

Distance  from  reference 
of  layer  at  node  J 

surface 

to  top 

(HJ1) 

31-40 

Wrap  angle 

(ANCLE1) 

*A11  Young's  moduli  must  be  scaled  down  by  a factor  of  10^ 

**The  last  card  of  the  set  must  contain  a number  greater  than 
the  total  node  number  (End  of  Data  Card). 

A pressurized  hemisphere-cylinder  shell  was  chosen  as  an  example 
for  data  preparation  for  the  computer  program.  The  shell  as  shown  in 
Figure  20  is  divided  into  four  segments.  There  are  two  boundary  points. 
With  the  boundary  condition  shown,  there  will  be  a boundary  release  at 
node  No.  1 in  the  z-direction  and  at  node  No.  50  in  the  r-direction. 

The  total  subtend  angles  for  segments  1 and  2 are  80°  and  10°  respec- 
tively (lines  7 and  9 of  Figure  21).  Note  that  the  angle  of  slope  between 
the  straight  line  segments  and  the  r-axis  for  segments  3 and  4 is  -90° 
since  the  node  numbers  are  increasing  downward  (see  line  13  of  Figure  21). 
A set  of  sample  data  cards  for  this  structure  is  shown  in  Figure  21. 
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Figure  20.  Hemisphere-Cylindrical  Shell. 
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